@                   Matching programmes for mediation analysis

Michael Lechner
Professor of Econometrics
SEW
University of St. Gallen
Switzerland

This programme estimates various direct and indirect effects under the assumption of no confounding on unobservable of treatment and mediators.

All rights reserved. Commercial use prohibited.

Using this programme for noncommercial purposes is free, if appropriate reference is given.

Version 1.0: Based on the matching procedures of bin_stat Version 3.1.1.

Oct., 8, 2013, M.Lechner

Decisions:
a) Common support: Is only based on extended model for propensity score (p(D=1|X,M) as it should lead to better predictions, i.e. more support problems than p(x)
b) Bootstrap and asymptotic inference for differences of effects: Based on i) Distribution of effects; ii) t-values computed under the assumption that effects are not correlated
c) Same tuning parameters used for both matching steps
d) Mediators can be included as additional variables in Mahalanobis distance


@

@ ========================================================================================================= @
@ ======================== BEGINN PROGRAMM ================================================================ @
@ ========================================================================================================= @

new;library pgraph;closeall;cls;ti1=date;

rndseed 12798306; outwidth 256;  

pfad    = "d:\\mlechner\\mediation1";	@ Path; there should exist be an /out and a /temp subdirecty;; check also path at end of programme             @ 

indat  = "data\\estimationnewspec";	    @ File that contains data input				        @

@ ========================================================================================================= @
@ ------------------ Controls for programme execution ---------------------------------------------------------------------- @
@ ========================================================================================================= @

step_1 	= 1;                @ final data set and descriptive stats of entry sample          @
step_3 	= 1;                @ Matching for all different combinations                       @
step_3a = 1;                @ .... corresponding descriptive statistics                     @
step_3b = 1;                @ ... score tests for probits                                   @
step_3c = 1;                @ descriptive statistics of the p-scores                        @
step_3d = 1;                @ statistics about common support; deleted and retained samples @
step_3e = 1;                @ statistics about match quality and distribution of weights    @
step_3ee = 1;               @ Correlation analysis of matching variables                    @

bootstrap = 999;            @ bootstrap = 0: no bootstrap; use 19, 99, 199, 499, 999 etc (!) as bootstrap sample sizes --> see MacKinnon 2006 @

number_boot = 999;         @ number of bootstrap replications for computing marginal effects in probit @
marg_eff_flag = 1;			@ if 0, marginal effects will not be computed @	

@ ========================================================================================================= @
@ ---------------- DEFINITIONs ------------------------------------ @
@ ========================================================================================================= @

outpfad  = pfad $+"\\out";

out_name = "HLM_BFL_6months";  	   					@ Output name for programm 	@

let treatment      = q24_23;					@ Treatmentvariable @
let treatment0     = treat0;   				@ 1-treatment; if it does not exist, generate in dataloop @

let mediators      = 
					//	basi3 pers3 spra3 info3 fach3     pvb3 
						basi6 pers6 spra6 info6 fach6     pvb6
					//	basi9 pers9 spra9 info9 fach9     pvb9
					//	basi12 pers12 spra12 info12 fach12 pvb12
;
// basisX personX sprachX infoX fachX pvbX with X=3,6,9,12 months.


let how_many_treatments = 1;    			@  treatment values higher than this will be deleted from the sample and a binary variable will be created         @

titel="Mediation";

let pers_nr = persnr;

let sample_w = 0 ;				@ sample weight; if this set to string 0, then the programme will add a sample weight with the name "WXXXX";@

cluster_boot = 1;							@ 1: Clustered bootstrap, jointly drawing variables with the same value of 'CLUSTER_VAR' @
let cluster_var = pbnrsiaw;


let outcomes =      	       
			arb1             arb2             arb3             arb4             arb5             arb6             arb7             arb8 
            arb9            arb10            arb11            arb12            arb13            arb14            arb15            arb16            arb17 
           arb18            arb19            arb20            arb21            arb22            arb23            arb24            arb25            arb26 
           arb27            arb28            arb29            arb30            arb31            arb32            arb33            arb34            arb35 
           arb36              
		   be1              be2              be3              be4              be5              be6             be7              be8              be9 
            be10             be11             be12             be13             be14             be15             be16             be17             be18 
            be19             be20             be21             be22             be23             be24             be25             be26             be27 
            be28             be29             be30             be31             be32             be33             be34             be35             be36 	
			stb_1 		stb_2 	stb_3	stb_4	stb_5	stb_6	stb_7	stb_8 	stb_9	stb_10
            stb_11		stb_12	stb_13	stb_14	stb_15  stb_16  stb_17  stb_18  stb_19  stb_20
            stb_21 		stb_22  stb_23  stb_24  stb_25  stb_26   stb_27   stb_28   stb_29   stb_30
            stb_31 		stb_32 	stb_33  stb_34  stb_35   stb_36 
			seb_1 		seb_2		seb_3		seb_4		seb_5	seb_6		seb_7	seb_8 	seb_9	seb_10
            seb_11	seb_12		seb_13		seb_14		seb_15  seb_16   	seb_17  seb_18   seb_19   seb_20
            seb_21  seb_22 	seb_23  	seb_24  	seb_25  seb_26   	seb_27   seb_28   seb_29   seb_30
            seb_31  seb_32 	seb_33  	seb_34  	seb_35  seb_36 
			atb1 atb2 		atb3		atb4		atb5	atb6		atb7	atb8 		atb9	atb10
            atb11	atb12	atb13		atb14		atb15  atb16   	atb17   atb18   atb19   atb20
            atb21  atb22  	atb23  	atb24  	atb25  atb26   	atb27   atb28   atb29   atb30
            atb31  atb32  	atb33  	atb34  	atb35  atb36 
			Stes1 	Stes2 	Stes3	Stes4	Stes5	Stes6	Stes7	Stes8 	Stes9	Stes10
            Stes11	Stes12	Stes13	Stes14	Stes15 Stes16   Stes17   Stes18   Stes19  Stes20
            Stes21 Stes22 Stes23 Stes24 Stes25 Stes26  Stes27   Stes28   Stes29   Stes30
            Stes31 Stes32 Stes33 Stes34 Stes35 Stes36 
			
 			mst_1 	mst_2 	mst_3	mst_4	mst_5	mst_6	mst_7	mst_8 	mst_9	mst_10
            mst_11	mst_12	mst_13	mst_14	mst_15   mst_16   mst_17   mst_18   mst_19   mst_20
            mst_21  mst_22  mst_23 mst_24  
		;
			
outcomes = mediators|outcomes;


let XVAR1    =
		frau      Q01      Q02      Q04      Q07    Q08_3    Q08_4    Q09_3    Q10_3    Q10_4    Q10_6    Q10_7    Q10_8    Q10_9   
		KANT_F   KANT_I   MUSP_R QUAL_UNG QUAL_ANG QUAL_OHN    ALQ_K   AL_2JS ET_SHARE    qmiss VERMIT_M VERMIT_L    
		F_Q01    F_Q02    F_Q04    F_Q07  F_Q08_3  F_Q08_4  F_Q09_3  F_Q10_3  F_Q10_4   
		F_FRAU F_MUSP_R F_QUALUN F_QUALAN F_QUALOH F_VMIT_L F_VMIT_M F_AL_2JS F_ET_SHA   F_ALQK    
		I_Q01    I_Q02    I_Q04    I_Q07  I_Q08_3  I_Q08_4  I_Q09_3  I_Q10_3  I_Q10_4 
		I_FRAU I_MUSP_R I_QUALUN I_QUALAN I_QUALOH I_VMIT_L I_VMIT_M I_AL_2JS I_ET_SHA   I_ALQK    /* std1     std2 */
		
		alters ALTER2S ZV_VERH bausland causland BCity MCity vverds ANZBE5S WZ_PRIM WZ_SEKU WZ_TERT
		FREMD_KS JB_SELBS JB_KADER JB_FACHF JB_HILFS GDP_PC
		;


XVAR2      = XVAR1|mediators;      				@ variables to be included in the specification tests of the propensity score specification @
                                      
let add_match_vars =frau;          		@ additional matching variables included in Mahalanobis distance ; they should be part of the control variables used in all probits @

add_match_vars_m = add_match_vars | mediators;

let score_weight = 5;			        	@ Values larger than one lead to an overweighting of the Score in the Mahalanobis-Metrik @			

match_check_vars = xvar1|xvar2|mediators;    @ variables used for descriptive stats about balancing @

x_extra_name 	= add_match_vars;           	@ variables used additionally to the score in the bias adjustment @
x_extra_name_M 	= add_match_vars_m;

let prozent_um_nachbar = 150;            	@ defines the radius in %-distance to 90%-quantil of distances obtained from one-to-one match @

let max_infl_weight = 5;                 	@ maximum share of weight of one observation compared to total weight; see HLW 2012; sampling weights ignored in this step; @

programm_test =0;         					@ if non-zero,  program uses X% random sample only @

let p_val_kernel = 1;						@ If = 0, empirical distribution function of t-statistic is used for bootstrap inference;
											  if = 1, p-values based on smoothed version as suggested by Racine and McKinnon (2007) @

let t_val_k = 1;							@ if = xxx; k-nearest neighbor estimation is used to compute the variance given the weights; k = xxxx * 2 * sqrt(N)	 
														Results do not appear sensitive to this choice
											  Instead of k-NN any type of nonparametric regression may be used, but the results are hardly affected by the choice of method at this stage@
		
	
@ ========================================================================================================= @
@ =================================== NO CHANGE BELOW THIS LINE OTHER THAN in DATALOOP ==================== @
@ ========================================================================================================= @
variablen_anzahl = 50;      @ Anzahl der Variablen, die bei Kreuztab in einem Schwung gelesen werdem    @
nr = 10000;					@ Rows per readr(.,nr)					                                    @
flag_on = 0;				@ 1: "output on";   0: "output reset" 	                                    @



let time_dimension_outcomes = 1;  @ define for how many periods exist; it is assumed that these the variables right to the first variables specified above @

com_var =   Xvar1;   test_var = xVar1|xVar2;
all_var = com_var|test_var;

outfile = 	out_name $+ ".out";
outdat1  =	out_name $+ "1",;   outdat2  = 	out_name $+ "2";    outdat3  = 	out_name $+ "3";    outdat4  = 	out_name $+ "4";    
outdat5  = 	out_name $+ "5";    outdat6  = 	out_name $+ "6";    outdat7  = 	out_name $+ "7";    outdat8  = 	out_name $+ "8";    outdat9  = 	out_name $+ "9";    

temp1="temptem1";temp2="temptem2";temp3="temptem3";temp4="temptem4";temp5="temptem5";
temp6="temptem6";temp7="temptem7"; temp8="temptem8";temp8_M="temptem8M";temp9="temptem9";temp10="temptem10";temp11="temptem11";
temp21="temptem21";temp22="temptem22";temp31="temptem31";temp41="temptem41";temp51="temptem51";temp61="temptem61";temp71="temptem71";temp81="temptem81";
temp21a="temptem21a";temp31a="temptem31a";temp41a="temptem41a";temp51a="temptem51a";temp61a="temptem61a";temp71a="temptem71a";temp81a="temptem81a";temp210="temp210";
temp9xx="tt9xx";temp99xx="tt99xx";

indat = pfad $+ "\\" $+ indat;  

outfile = outpfad $+ "\\" $+ outfile; outdat1  = outpfad $+ "\\" $+ outdat1;outdat2  = outpfad $+ "\\" $+ outdat2;outdat3  = outpfad $+ "\\" $+ outdat3;
outdat4  = outpfad $+ "\\" $+ outdat4;outdat5  = outpfad $+ "\\" $+ outdat5;outdat6  = outpfad $+ "\\" $+ outdat6;outdat7  = outpfad $+ "\\" $+ outdat7;
outdat8  = outpfad $+ "\\" $+ outdat8;outdat9  = outpfad $+ "\\" $+ outdat9;
temppfad = pfad $+"\\temp";

temp1=temppfad $+ "\\" $+temp1;temp2=temppfad $+ "\\" $+temp2;temp3=temppfad $+ "\\" $+temp3;
temp4=temppfad $+ "\\" $+temp4;temp5=temppfad $+ "\\" $+temp5;temp6=temppfad $+ "\\" $+temp6;
temp7=temppfad $+ "\\" $+temp7;temp8=temppfad $+ "\\" $+temp8;temp8_M=temppfad $+ "\\" $+temp8_M;temp9=temppfad $+ "\\" $+ temp9;
temp10=temppfad $+ "\\" $+ temp10;temp11=temppfad $+ "\\" $+ temp11;
temp21=temppfad $+ "\\" $+ temp21;temp22=temppfad $+ "\\" $+ temp22;temp31=temppfad $+ "\\" $+ temp31;temp41=temppfad $+ "\\" $+ temp41;temp51=temppfad $+ "\\" $+ temp51;
temp61=temppfad $+ "\\" $+ temp61;temp71=temppfad $+ "\\" $+ temp71;temp81=temppfad $+ "\\" $+ temp81;
temp21a=temppfad $+ "\\" $+ temp21a;temp31a=temppfad $+ "\\" $+ temp31a;temp41a=temppfad $+ "\\" $+ temp41a;temp51a=temppfad $+ "\\" $+ temp51a;
temp61a=temppfad $+ "\\" $+ temp61a;temp71a=temppfad $+ "\\" $+ temp71a;temp81a=temppfad $+ "\\" $+ temp81a;temp210=temppfad $+ "\\" $+ temp210;
temp9xx=temppfad $+ "\\" $+ temp9xx;temp99xx=temppfad $+ "\\" $+ temp99xx;

load path=^outpfad;	save path=^outpfad;    
if flag_on == 1; output file=^outfile on; else; output file=^outfile reset;  endif;

if bootstrap < 0; bootstrap = 0;  endif;

@ ========================================================================================================= @
@ ------ 1. Descriptive statistics of input files --------------------------------------------------------- @
@ ========================================================================================================= @



if step_1 == 1;
 format 16,16; 
	"******************************** Dateninput ********************************";
	call dstatneu(indat);   
 
	dataloop ^indat ^temp22;  /* this dataloop has to be modified in an application */

		make basin1 = basisn1;              
		make persn1 = personn1;             
		make spran1 = sprachn1;             

        make basi3 = basis3;
		make pers3 = person3;
		make spra3 = sprach3;
		
		make basi6 = basis6;
		make pers6 = person6;
		make spra6 = sprach6;

        make basi9 = basis9;
		make pers9 = person9;
		make spra9 = sprach9;
		
		
		make basi12 = basis12;
		make pers12 = person12;
		make spra12 = sprach12;
		
		make std1   = startdate == 1;
		make std2   = startdate == 2;
		make treat0 = 1 - q24_23;
 
		drop
			be_1, 	be_2 ,		be_3,		be_4,		be_5,	be_6	,	be_7	,be_8 ,	be_9,	be_10,
            be_11,	be_12,		be_13,		be_14,		be_15 , be_16 ,  	be_17 , be_18  , be_19 ,  be_20,
            be_21,  be_22 , 	be_23,  	be_24,  	be_25 , be_26  , 	be_27  , be_28 ,  be_29 ,  be_30,
            be_31,  be_32 , 	be_33 , 	be_34 , 	be_35,  be_36 ,  	be_37,   be_38 ,  be_39 ,  be_40  ,          
			be_41,  be_42,  	be_43 , 	be_44 , 	be_45 , be_46   ,	be_47 ,
			arb_1, 	arb_2 ,		arb_3,		arb_4,		arb_5,	arb_6	,	arb_7	,arb_8 ,	arb_9,	arb_10,
            arb_11,	arb_12,		arb_13,		arb_14,		arb_15 , arb_16 ,  	arb_17 , arb_18  , arb_19 ,  arb_20,
            arb_21,  arb_22 , 	arb_23,  	arb_24,  	arb_25 , arb_26  , 	arb_27  , arb_28 ,  arb_29 ,  arb_30,
            arb_31,  arb_32 , 	arb_33 , 	arb_34 , 	arb_35,  arb_36 ,  	arb_37,   arb_38 ,  arb_39 ,  arb_40  ,          
			arb_41,  arb_42,  	arb_43 , 	arb_44 , 	arb_45 , arb_46   ,	arb_47 ;
 
 		drop	
			  	seb_37,   seb_38 ,  seb_39 ,  seb_40  ,          
			seb_41,  seb_42,  	seb_43 , 	seb_44 , 	seb_45 , seb_46   ,	seb_47 ,
			atb_1 ,	atb_2 ,		atb_3,		atb_4,		atb_5,	atb_6	,	atb_7	,atb_8 	,	atb_9,	atb_10,
            atb_11,	atb_12,		atb_13,		atb_14,		atb_15 , atb_16  , 	atb_17   ,atb_18 ,  atb_19  , atb_20,
            atb_21,  atb_22 , 	atb_23,  	atb_24 , 	atb_25 , atb_26 ,  	atb_27  , atb_28 ,  atb_29,   atb_30,
            atb_31,  atb_32,  	atb_33 , 	atb_34 , 	atb_35  ,atb_36 ,  	atb_37,   atb_38 ,  atb_39 ,  atb_40  ,          
			atb_41,  atb_42,  	atb_43   	atb_44 , 	atb_45  ,atb_46  , 	atb_47;
		drop Stes_1 	Stes_2 	Stes_3	Stes_4	Stes_5	Stes_6	Stes_7	Stes_8 	Stes_9	Stes_10
	           Stes_11	Stes_12	Stes_13	Stes_14	Stes_15 Stes_16   Stes_17   Stes_18   Stes_19   Stes_20
            Stes_21 Stes_22 Stes_23 Stes_24 Stes_25 Stes_26  Stes_27   Stes_28   Stes_29   Stes_30
            Stes_31 Stes_32 Stes_33 Stes_34 Stes_35 Stes_36 
		mst_25   mst_26   mst_27   mst_28   mst_29   mst_30
            mst_31 mst_32 mst_33 mst_34 mst_35 mst_36 	;
		drop
			  Stes_37  , Stes_38 ,  Stes_39  , Stes_40  ,          
			Stes_41, Stes_42 ,Stes_43 ,Stes_44 ,Stes_45 ,Stes_46  , Stes_47,
			lzt_1 ,	lzt_2 ,	lzt_3,	lzt_4,	lzt_5,	lzt_6	,lzt_7,	lzt_8 ,	lzt_9	,lzt_10,
            lzt_11,	lzt_12,	lzt_13,	lzt_14,	lzt_15 , lzt_16 ,  lzt_17 ,  lzt_18  , lzt_19  , lzt_20,
            lzt_21 , lzt_22 , lzt_23 , lzt_24 , lzt_25  ,lzt_26,   lzt_27  , lzt_28 ,  lzt_29 ,  lzt_30,
            lzt_31 , lzt_32 , lzt_33 ,	lzt_34 , lzt_35  ,lzt_36 ,  lzt_37 ,  lzt_38 ,  lzt_39  , lzt_40 ,           
			lzt_41 , lzt_42 , lzt_43  ,lzt_44 , lzt_45 , lzt_46 ,  lzt_47;
		drop	
			aus_1 ,	aus_2 ,	aus_3,	aus_4,	aus_5,	aus_6,	aus_7	,aus_8 	,aus_9	,aus_10,
            aus_11,	aus_12,	aus_13,	aus_14,	aus_15  ,aus_16 ,  aus_17  , aus_18  , aus_19 ,  aus_20,
            aus_21 , aus_22,  aus_23 , aus_24  ,aus_25  ,aus_26 ,  aus_27 ,  aus_28  , aus_29  , aus_30,
            aus_31,  aus_32 , aus_33 , aus_34 , aus_35 , aus_36  , aus_37 ,  aus_38 ,  aus_39 ,  aus_40 ,           
			aus_41 , aus_42,  aus_43 , aus_44 , aus_45 , aus_46 ,  aus_47;
	drop
			  mst_37  , mst_38 ,  mst_39  , mst_40 ,           
			mst_41 , mst_42 , mst_43 , mst_44   mst_45  , mst_46  , mst_47,
			  stb_37 ,  stb_38  , stb_39  , stb_40,            
			stb_41 , stb_42 , stb_43 , stb_44 , stb_45 ,  stb_46 ,  stb_47	;			
  
		drop 	be37,   be38 ,  be39 ,  be40  ,          
			be41,  be42,  	be43 , 	be44 , 	be45 , be46   ,	be47 ,
			arb37,   arb38 ,  arb39 ,  arb40  ,          
			arb41,  arb42,  	arb43 , 	arb44 , 	arb45 , arb46   ,	arb47 ;
 
 		drop	
			seb1, 	seb2 ,		seb3,		seb4,		seb5,	seb6	,	seb7	,seb8 ,	seb9,	seb10,
            seb11,	seb12,		seb13,		seb14,		seb15 , seb16 ,  	seb17 , seb18  , seb19 ,  seb20,
            seb21,  seb22 , 	seb23,  	seb24,  	seb25 , seb26  , 	seb27  , seb28 ,  seb29 ,  seb30,
            seb31,  seb32 , 	seb33 , 	seb34 , 	seb35,  seb36 ,  	seb37,   seb38 ,  seb39 ,  seb40  ,          
			seb41,  seb42,  	seb43 , 	seb44 , 	seb45 , seb46   ,	seb47 ,
			atb37,   atb38 ,  atb39 ,  atb40  ,          
			atb41,  atb42,  	atb43   	atb44 , 	atb45  ,atb46  , 	atb47;
		drop	
			Stes37  , Stes38 ,  Stes39  , Stes40  ,          
			Stes41, Stes42 ,Stes43 ,Stes44 ,Stes45 ,Stes46  , Stes47,
			lzt1 ,	lzt2 ,	lzt3,	lzt4,	lzt5,	lzt6	,lzt7,	lzt8 ,	lzt9	,lzt10,
            lzt11,	lzt12,	lzt13,	lzt14,	lzt15 , lzt16 ,  lzt17 ,  lzt18  , lzt19  , lzt20,
            lzt21 , lzt22 , lzt23 , lzt24 , lzt25  ,lzt26,   lzt27  , lzt28 ,  lzt29 ,  lzt30,
            lzt31 , lzt32 , lzt33 ,	lzt34 , lzt35  ,lzt36 ,  lzt37 ,  lzt38 ,  lzt39  , lzt40 ,           
			lzt41 , lzt42 , lzt43  ,lzt44 , lzt45 , lzt46 ,  lzt47;
		drop	
			aus1 ,	aus2 ,	aus3,	aus4,	aus5,	aus6,	aus7	,aus8 	,aus9	,aus10,
            aus11,	aus12,	aus13,	aus14,	aus15  ,aus16 ,  aus17  , aus18  , aus19 ,  aus20,
            aus21 , aus22,  aus23 , aus24  ,aus25  ,aus26 ,  aus27 ,  aus28  , aus29  , aus30,
            aus31,  aus32 , aus33 , aus34 , aus35 , aus36  , aus37 ,  aus38 ,  aus39 ,  aus40 ,           
			aus41 , aus42,  aus43 , aus44 , aus45 , aus46 ,  aus47;
	drop
			mst1 ,	mst2 ,	mst3,	mst4,	mst5,	mst6,	mst7,	mst8 ,	mst9	,mst10,
            mst11,	mst12,	mst13,	mst14	mst15 ,  mst16 ,  mst17 ,  mst18 ,  mst19,   mst20,
            mst21 , mst22,  mst23 , mst24 , mst25  , mst26  , mst27  , mst28 ,  mst29 ,  mst30,
            mst31,  mst32 , mst33 , mst34 , mst35 ,  mst36 ,  mst37  , mst38 ,  mst39  , mst40 ,           
			mst41 , mst42 , mst43 , mst44   mst45  , mst46  , mst47,
			stb1 ,	stb2 ,	stb3,	stb4,	stb5	,stb6	stb7,	stb8 ,	stb9,	stb10,
            stb11,	stb12,	stb13,	stb14,	stb15  , stb16 ,  stb17  , stb18  , stb19  , stb20,
            stb21 , stb22,  stb23 , stb24 , stb25 ,  stb26 ,  stb27 ,  stb28 ,  stb29  , stb30,
            stb31 , stb32 , stb33 , stb34 , stb35 ,  stb36 ,  stb37 ,  stb38  , stb39  , stb40,            
			stb41 , stb42 , stb43 , stb44 , stb45 ,  stb46 ,  stb47	;	
	endata;     
	
	if sample_w $== 0;  	sample_w = add_sample_weight(temp22,temp1,nr); /* adds weight to TEMP22; temp1 is only temporary */  
							call dstatneu(temp22);
	else;				 	call dstatneu(temp22); call dstatneu_w(temp22,getname_neu(temp22),sample_w);   	  				 
	endif;
endif;

/*
x=readvar(temp22,outcomes);
m=meanc(x);
s=stdc(x);
nullx= (s .< 1e-8) .and (m .< 1e-8);

for i(1,rows(outcomes),1);
	format 8,8; $outcomes[i];;format 16,8;m[i];;s[i];; format 16,0;nullx[i];
	
endfor;	
end;
*/


temp21alt = temp21; 
No_w = stdc(readvar(temp22,sample_w)) == 0;		/* contains no variation in weights */

for b_i(0,bootstrap,1);
	temp21 = temp21alt;
	if b_i == 1;"iteration  ";; endif;
	output on; format 1,0; b_i;;
        
	if b_i > 0; 
		if cluster_boot == 1;
			call exctsmpl_cluster_100(temp22,temp210,cluster_var,10000,1000,unique(readvar(temp22,cluster_var),1));
		else;				
			call exctsmpl(temp22,temp210,100);                  /* zweites file ist output */
		endif;
		call unique_pers_nummer(pers_nr,temp210,temp21);  /* zweites file ist output */
	
	else;                     temp21 = temp22;
	endif;       
	temp1 = temp21; 
    
  @ ========================================================================================================= @
  @ ------ 3. Matching for all pairwise comparisons (NOT TESTED FOR MULTIPLE TREATMENTS --------------------------------------------------------- @
  @ ========================================================================================================= @
  
	if step_3 == 1;
		for i(1,how_many_treatments,1);  
			if b_i == 0; " ****************************** Reference treatment ";; format 1,0; i;; " *****************************************";endif;
			for j(0,how_many_treatments,1);
				if i /= j;
					addvar = -1;
					if addvar == -1;  indep_name = com_var;  if b_i == 0;   "No additional variables specified"; endif; 
					else;             indep_name = com_var | addvar; 
					endif;
   
					indep_name_M = indep_name | mediators;
   
					{dep_name} = get_data_for_pair(temp1,temp2,treatment,i,j,nr);  @ selected sample in temp2 @
					if b_i == 0;
						?;"===================================================================================================";
						"  Reference population ";; format 3,0; i;; " *** Control population: ";; j;
						"---------------------------------------------------------------------------------------------------"; ?;
					endif;
					if (step_3a == 1) and (b_i == 0); 
						call dstatneu(temp2); 
						call new_kreuztab(temp2,treatment,getname_neu(temp2),nr,variablen_anzahl); 
					endif;
					if b_i == 0; ausgabepr = 1; call new_kreuztab(temp2,treatment,dep_name|indep_name,nr,variablen_anzahl); else; ausgabepr = 0; endif;
					if no_w;{b,ll,df,p,m,hxss,opg} 					= probit_a(temp2,dep_name,indep_name,  zeros(rows(indep_name)+1,1),1,1e-8,1,nr,ausgabepr);
							{b_m,ll_m,df_m,p_m,m_m,hxss_m,opg_m} 	= probit_a(temp2,dep_name,indep_name_m,zeros(rows(indep_name_m)+1,1),1,1e-8,1,nr,ausgabepr);
							
					else;	{b,ll,df,p,m,hxss,opg} 					= probit_a_weight(temp2,dep_name,indep_name,  zeros(rows(indep_name)+1,1),1,1e-8,1,nr,ausgabepr,sample_w);
							{b_m,ll_m,df_m,p_m,m_m,hxss_m,opg_m} 	= probit_a_weight(temp2,dep_name,indep_name_m,zeros(rows(indep_name_m)+1,1),1,1e-8,1,nr,ausgabepr,sample_w);
					endif;
					if b_i == 0;
						b_probit_0 = b;
						p_probit_0 = p;
						b_probit_0_m = b_m;
						p_probit_0_m = p_m;
						if marg_eff_flag;						
							call marg_effekte(temp2,temp9xx,temp99xx,dep_name,indep_name,  sample_w,2,number_boot,cluster_var,cluster_boot);
							?;"Including mediators";
							call marg_effekte(temp2,temp9xx,temp99xx,dep_name,indep_name_M,sample_w,2,number_boot,cluster_var,cluster_boot);
						endif;				
					else;
						if b[rows(b)] == 0;			/* probit collapsed */
							b = b_probit_0;
							p = p_probit_0;	
							b_m = b_probit_0_m;
							p_m = p_probit_0_m;	
						endif;
				    endif;	
					pt 		= p;
					pt_m 	= p_m;
					nix=readvar(temp2,dep_name); if minc(nix)/=0 or maxc(nix)/=1; "fdehler"; end; endif;
					if b_i == 0;
						_ptek = outpfad $+"\\pscore_" $+ titel $+ "_D0.tkf"; call histp(selif(p,nix.==0),100); 
						_ptek = outpfad $+"\\pscore_" $+ titel $+ "_D1.tkf"; call histp(selif(p,nix.==1),100); 
						_ptek = outpfad $+"\\pscore_" $+ titel $+ "_M_D0.tkf"; call histp(selif(p_m,nix.==0),100); 
						_ptek = outpfad $+"\\pscore_" $+ titel $+ "_M_D1.tkf"; call histp(selif(p_m,nix.==1),100); 
					endif;
					if (step_3b == 1) and (b_i == 0) and no_w;
						call omitted_bin(b,hxss,m,opg,temp2,dep_name,indep_name,temp2,xVar2,nr,2,0,0,0,0,0,0,0,0,0,0); 
						call Jar_Bera(temp2,dep_name,indep_name,b,nr); 
						call heteros_bin(b,hxss,m,opg,temp2,dep_name,indep_name,temp2,indep_name,nr,2,0,0,0,0,0,0,0,0,0,0);
						?; " Including mediators";	call omitted_bin(b_m,hxss_m,m_m,opg_m,temp2,dep_name,indep_name_m,temp2,xVar2,nr,2,0,0,0,0,0,0,0,0,0,0); 
						?; " Including mediators";  call Jar_Bera(temp2,dep_name,indep_name_m,b_m,nr); 
						?; " Including mediators";	call heteros_bin(b_m,hxss_m,m_m,opg_m,temp2,dep_name,indep_name_m,temp2,indep_name_m,nr,2,0,0,0,0,0,0,0,0,0,0);
					endif;                
					let p_score = xb ;
					let p_score_M = xb_M;
					{score_name_all_t}   = add_xb_p_multiple(temp1,temp3,temp9,p_score,indep_name,b,i,j,nr);
					{score_name_all_t_M} = add_xb_p_multiple(temp3,temp3,temp9,p_score_M,indep_name_M,b_M,i,j,nr);
					
					if (j == 0) or ((j == 1) and (i == 0)); score_name_all 		=                  score_name_all_t;
															score_name_all_M 	=                  score_name_all_t_M;
					else;                                   score_name_all 		= score_name_all | score_name_all_t;
															score_name_all_M 	= score_name_all_M | score_name_all_t_M;
					endif;
				endif; /* i /= j; */
			endfor;  @ j TEMP3 contains the relevant data plus the scores !@
 
			if (step_3c == 1) and (b_i == 0);     call dstatneu(temp3);    call new_kreuztab(temp3,treatment,getname_neu(temp3),nr,variablen_anzahl); endif;
			{share_remaining} = commonsupport_max_treated_m_both(temp3,temp4,temp8,temp9,treatment,score_name_all_M,i,nr,max_infl_weight); @ temp 4 enth�lt gute daten; temp5 schlechte daten @
			/* common support is based on the larger model as it should predict better */
		
			?;?;"***************************************************************************";
			"------------- remaining in % -----------------";; format 4,0; (share_remaining*100);
			if (step_3d == 1) and (b_i == 0); 
				"------------- deskr. stats of sample to be continued with --------------";
				"------------- remaining observations -----------------";
				if no_w; call dstatneu(temp4);        call new_kreuztab(temp4,treatment,treatment|indep_name,nr,variablen_anzahl);
				else;	 call dstatneu_w(temp4,getname_neu(temp4),sample_w);
						 call new_kreuztab_w(temp4,treatment,treatment|indep_name,sample_w);
				endif;
		 
				"------------- deleted  sample -----------------";
				if share_remaining < 1;
					if no_w; 	call dstatneu(temp9);							call new_kreuztab(temp9,  treatment,treatment|indep_name,nr,variablen_anzahl);
					else;		call dstatneu_w(temp9,getname_neu(temp9),sample_w);	call new_kreuztab_w(temp9,treatment,treatment|indep_name,sample_w);
					endif;
				else; "No observations deleted";
				endif;
			endif;
			
			@ ========================================================================================================= @
			if b_i == 0;" ***********************  NOW Matching on the common support *******";endif;
			@ ========================================================================================================= @
			for j(0,how_many_treatments,1);
				if i /= j;
					addvar = -1;
					addvar2 = -1;
					@ Here comes loop to get to the ATE by estimated ATENT and ATET and weighting them by pop.schares of instrument @
					for ate_i(0,1,1);
						if     ate_i == 0;         instru = treatment0;  /* instrument umgedreht um -ATNT zu sch�tzen */
						elseif ate_i == 1;         instru = treatment;   
						endif;
						{dep_name} = get_data_for_pair(temp4,temp5,instru,i,j,nr);  @ selected sample in temp2 @
						if b_i == 0;
							?;"===================================================================================================";
							"  Reference population ";; format 3,0; i;; " *** Control population: ";; j;
							"---------------------------------------------------------------------------------------------------"; ?;
						endif;
						@ Radius Matching with 2 steps @
						if addvar2 == -1;	atc = add_match_vars; 			atc_m = add_match_vars_m; 
						else; 				atc = add_match_vars|addvar2; 	atc_m = add_match_vars_m|addvar2; 
						endif;     
						if i < j; 	p_score_name_t = score_name_all[j]; 	p_score_name_t_m = score_name_all_m[j]; 
						else; 		p_score_name_t = score_name_all[j+1];   p_score_name_t_m = score_name_all_m[j+1];
						endif;
						if b_i == 0; "Propensity scores and variables";; 			format 8,8;$p_score_name_t;; $atc'; endif;
						if b_i == 0; "Propensity scores and variables Mediators";; 	format 8,8;$p_score_name_t_M;; $atc_M'; endif;
						if (step_3ee == 1) and (b_i == 0) and no_w; 
											"========== P(X) ======================";
							if atc == 0; 	call correlation_analysis_matching(temp5,p_score_name_t,dep_name);
							else;		 	call correlation_analysis_matching(temp5,p_score_name_t_m|atc_m,dep_name);
							endif;		
							if atc_m == 0;	"========== P(X,M) ======================";
											call correlation_analysis_matching(temp5,p_score_name_t_m,dep_name);
							else;		 	call correlation_analysis_matching(temp5,p_score_name_t_m|atc_m,dep_name);
							endif;		
						endif;
						if no_w;	threadstat {pn1,  pn0,  w_mit_nuller}  = match_p_and_time_constant_m_s(  temp5,p_score_name_t,  dep_name,pers_nr,atc,	score_weight,prozent_um_nachbar);
									threadstat {pn1_m,pn0_m,w_mit_nuller_m}= match_p_and_time_constant_m_s(  temp5,p_score_name_t_m,dep_name,pers_nr,atc_m,score_weight,prozent_um_nachbar);
									threadjoin;
						else;		threadstat {pn1,  pn0,  w_mit_nuller}  = match_p_and_time_constant_m_s_w(temp5,p_score_name_t,  dep_name,pers_nr,atc,  score_weight,prozent_um_nachbar,sample_w);
									threadstat {pn1_M,pn0_M,w_mit_nuller_M}= match_p_and_time_constant_m_s_w(temp5,p_score_name_t_m,dep_name,pers_nr,atc_m,score_weight,prozent_um_nachbar,sample_w);
									threadjoin;
						endif;
						"Scoreweight: ";; format 8,3; score_weight;  @ pn0 enthaelt auch die Gewichte, nicht normiert @
						if (step_3e == 1) and (b_i == 0);
							if addvar == -1; 	mcv 	= match_check_vars; 
												mcv_m 	= match_check_vars; 
							else;  				mcv 	= match_check_vars  |addvar; 
												mcv_M 	= match_check_vars|addvar; 
							endif;   
						endif;
						y = readvar(temp5,dep_name); if b_i == 0;"Hier muss 0 rauskommen, wenn die Ordnung ok ist: ";; format 8,4; sumc( y .* w_mit_nuller);; sumc( y .* w_mit_nuller_m); endif;					
						 						
						if no_w;		w_komplett   = w_mit_nuller   +   y/sumc(y);  /* setzt die gewichte der treated auf 1/n1 */ 
										w_komplett_M = w_mit_nuller_m + y/sumc(y); 
						else;			w = readvar(temp5,sample_w); 
										w_komplett   = w_mit_nuller   + y .* w;
										w_komplett_M = w_mit_nuller_m + y .* w;
										w_komplett 		= y .* w_komplett 	./ sumc(y .* w_komplett) 	+ (1-y) .* w_komplett 	./ sumc((1-y) .* w_komplett);
										w_komplett_m 	= y .* w_komplett_m ./ sumc(y .* w_komplett_m) 	+ (1-y) .* w_komplett_m ./ sumc((1-y) .* w_komplett_m);
						endif;		
						{w,treated,w_name} 		= get_pn_into_w_m_s(temp8,temp5,w_komplett,pers_nr,dep_name,nr);
						{w_m,treated,w_name_M} 	= get_pn_into_w_m_s(temp8_M,temp5,w_komplett_M,pers_nr,dep_name,nr);
						
						if (step_3e == 1) and (b_i == 0);
							call w_auswertung_simple(w,  treated);		call w_auswertung_simple(w_M,treated);
						endif;
						@ ========================================================================================================= @
						names = getname_neu(temp8);						
						
						{nix,o_i}= indices(temp8,outcomes);
						if time_dimension_outcomes > 1;
							for jj(1,rows(outcomes),1);
								index_j = o_i[jj];
								all_outcomes_t = names[index_j:(index_j+time_dimension_outcomes-1)];
								if jj == 1; all_outcomes = all_outcomes_t; else; all_outcomes = all_outcomes | all_outcomes_t;  endif;
							endfor;       
						else;	
							all_outcomes = outcomes;
						endif;	

						if (((i == 0) and (j == 1)) or (j == 0)); create_yes = 1; else; create_yes = 0;  endif;

						if ate_i == 0;  /* ------ ATENT --------------- */ /* note: fuer ate_i == 0 sind die rolen der gruppe 0 und 1 vertauscht */                           
							outdat_t1 = outdat1 $+ "_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */
							outdat_t2 = outdat2 $+ "_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */
							outdat_t3 = outdat3 $+ "_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */

							outdat_t1_M = outdat1 $+ "_M_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */
							outdat_t2_M = outdat2 $+ "_M_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */
							outdat_t3_M = outdat3 $+ "_M_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */

							if no_w;
								threadstat {y1xt, y0xt, y1xt_s, y0xt_s, y0xt_975, y1xt_975, y0xt_025, y1xt_025, n1x,n0x, n_refx,w_match_WOLSN}  	=evaluate_m_s_all_2(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,1,t_val_k); 
								threadstat {y1xt2,y0xt2,y1xt_s2,y0xt_s2,y0xt_9752,y1xt_9752,y0xt_0252,y1xt_0252,n1x2,n0x2,n_refx2,w_match2N}		=evaluate_m_s_all_2(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,2,t_val_k); 
								threadstat {y1xt3,y0xt3,y1xt_s3,y0xt_s3,y0xt_9753,y1xt_9753,y0xt_0253,y1xt_0253,n1x3,n0x3,n_refx3,w_match3N}		=evaluate_m_s_all_2(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,0,t_val_k); 
								threadjoin; 

								threadstat {y1xt_M, y0xt_M, y1xt_s_M, y0xt_s_M, y0xt_975_M, y1xt_975_M, y0xt_025_M, y1xt_025_M, n1x_M, n0x_M, n_refx_M, w_match_WOLSN_M}=evaluate_m_s_all_2(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,1,t_val_k); 
								threadstat {y1xt2_M,y0xt2_M,y1xt_s2_M,y0xt_s2_M,y0xt_9752_M,y1xt_9752_M,y0xt_0252_M,y1xt_0252_M,n1x2_M,n0x2_M,n_refx2_M,w_match2N_M}	=evaluate_m_s_all_2(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,2,t_val_k); 
								threadstat {y1xt3_M,y0xt3_M,y1xt_s3_M,y0xt_s3_M,y0xt_9753_M,y1xt_9753_M,y0xt_0253_M,y1xt_0253_M,n1x3_M,n0x3_M,n_refx3_M,w_match3N_M}	=evaluate_m_s_all_2(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,0,t_val_k); 
								threadjoin;

							else;
								threadstat {y1xt, y0xt, y1xt_s, y0xt_s, y0xt_975, y1xt_975, y0xt_025, y1xt_025, n1x,n0x, n_refx,w_match_WOLSN,w1N}  =evaluate_m_s_all_2w(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,1,sample_w,t_val_k); 
								threadstat {y1xt2,y0xt2,y1xt_s2,y0xt_s2,y0xt_9752,y1xt_9752,y0xt_0252,y1xt_0252,n1x2,n0x2,n_refx2,w_match2N,w12N}	=evaluate_m_s_all_2w(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,2,sample_w,t_val_k); 
								threadstat {y1xt3,y0xt3,y1xt_s3,y0xt_s3,y0xt_9753,y1xt_9753,y0xt_0253,y1xt_0253,n1x3,n0x3,n_refx3,w_match3N,w13N}	=evaluate_m_s_all_2w(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,0,sample_w,t_val_k); 
								threadjoin;						
								
								threadstat {y1xt_M, y0xt_M, y1xt_s_M, y0xt_s_M, y0xt_975_M, y1xt_975_M, y0xt_025_M, y1xt_025_M, n1x_M,n0x, n_refx_M,w_match_WOLSN_M,w1N_M}  =evaluate_m_s_all_2w(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,1,sample_w,t_val_k); 
								threadstat {y1xt2_M,y0xt2_M,y1xt_s2_M,y0xt_s2_M,y0xt_9752_M,y1xt_9752_M,y0xt_0252_M,y1xt_0252_M,n1x2_M,n0x2,n_refx2_M,w_match2N_M,w12N_M}	=evaluate_m_s_all_2w(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,2,sample_w,t_val_k); 
								threadstat {y1xt3_M,y0xt3_M,y1xt_s3_M,y0xt_s3_M,y0xt_9753_M,y1xt_9753_M,y0xt_0253_M,y1xt_0253_M,n1x3_M,n0x3,n_refx3_M,w_match3N_M,w13N_M}	=evaluate_m_s_all_2w(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,0,sample_w,t_val_k); 
								threadjoin;						
							endif;
							if b_i == 0; 
								?;output on;
								"========================================================================";
								"====== Balancing Test before matching ATENT xx====";
								if no_w;	?;"P(X)";	call test_quality_of_matched_s_m_s(temp8,mcv|score_name_all,dep_Name,1); 
											?;"P(X,M)";	call test_quality_of_matched_s_m_s(temp8_m,mcv_m|score_name_all_m,dep_Name,1); 
								else;		?;"P(X)";   w0x = readvar(temp8,sample_w);
								                        dxx = readvar(temp8,dep_name);
											            call test_quality_of_matched_s_m_s_w(temp8,mcv|score_name_all,dep_Name,selif(w0x,dxx .== 0),w1N); 
											?;"P(X,M)"; w0x_m = readvar(temp8_m,sample_w);
											            dxx_M = readvar(temp8_m,dep_name);
											call test_quality_of_matched_s_m_s_w(temp8_m,mcv_m|score_name_all_m,dep_Name,selif(w0x_m,dxx_M .== 0),w1N_m); 
								endif;	
								"========================================================================";
								"====== Balancing Test with weights from radius matching ATENT ====";
								if no_w;	?;"P(X)";	call test_quality_of_matched_s_m_s(temp8,mcv|score_name_all,dep_Name,w_match3N); 
											?;"P(X,M)";	call test_quality_of_matched_s_m_s(temp8_m,mcv_m|score_name_all_m,dep_Name,w_match3N_m); 
								else;		?;"P(X)";	call test_quality_of_matched_s_m_s_w(temp8,mcv|score_name_all,dep_Name,w_match3N,w13N); 
											?;"P(X,M)";	call test_quality_of_matched_s_m_s_w(temp8_m,mcv_m|score_name_all_m,dep_Name,w_match3N_m,w13N_m);
								endif;
								"========================================================================";
								"====== Balancing Test with weights from linear bias adjustment ATENT ====";
								if no_w;	?;"P(x)";	call test_quality_of_matched_s_m_s(temp8,mcv|score_name_all,dep_Name,w_match_WOLSN); 
											?;"P(x,M)";	call test_quality_of_matched_s_m_s(temp8_m,mcv_m|score_name_all_m,dep_Name,w_match_WOLSN_m); 
								else;		?;"P(X)";	call test_quality_of_matched_s_m_s_w(temp8,mcv|score_name_all,dep_Name,w_match_WOLSN,w1N);
											?;"P(X,M)";	call test_quality_of_matched_s_m_s_w(temp8,mcv_m|score_name_all_m,dep_Name,w_match_WOLSN_m,w1N_m); 	
								endif;		
								"========================================================================";
								call write_outcomes_in_file(outdat_t1,y0xt, y1xt, y0xt_s, y1xt_s, y1xt_975, y0xt_975, y1xt_025, y0xt_025, n0x, n1x, n_refx,share_remaining,all_outcomes,create_yes); 
								call write_outcomes_in_file(outdat_t2,y0xt2,y1xt2,y0xt_s2,y1xt_s2,y1xt_9752,y0xt_9752,y1xt_0252,y0xt_0252,n0x2,n1x2,n_refx2,share_remaining,all_outcomes,create_yes); 
								call write_outcomes_in_file(outdat_t3,y0xt3,y1xt3,y0xt_s3,y1xt_s3,y1xt_9753,y0xt_9753,y1xt_0253,y0xt_0253,n0x3,n1x3,n_refx3,share_remaining,all_outcomes,create_yes);
				
								call write_outcomes_in_file(outdat_t1_M,y0xt_m, y1xt_m, y0xt_s_m, y1xt_s_m, y1xt_975_m, y0xt_975_m, y1xt_025_m, y0xt_025_m, n0x_m, n1x_m, n_refx_m,share_remaining,all_outcomes,create_yes); 
								call write_outcomes_in_file(outdat_t2_M,y0xt2_m,y1xt2_m,y0xt_s2_m,y1xt_s2_m,y1xt_9752_m,y0xt_9752_m,y1xt_0252_m,y0xt_0252_m,n0x2_m,n1x2_m,n_refx2_m,share_remaining,all_outcomes,create_yes); 
								call write_outcomes_in_file(outdat_t3_m,y0xt3_m,y1xt3_m,y0xt_s3_m,y1xt_s3_m,y1xt_9753_m,y0xt_9753_m,y1xt_0253_m,y0xt_0253_m,n0x3_m,n1x3_m,n_refx3_m,share_remaining,all_outcomes,create_yes);
				
							endif;
							if (b_i == 0) or (b_i == 1); 
								threadbegin;allxthetas  =              (y1xt  - y0xt)';     allxstds  =              sqrt((y1xt_s .^2) + (y0xt_s  .^2))'; threadend; 
								threadbegin;allxthetas2 =              (y1xt2 - y0xt2)';    allxstds2 =              sqrt((y1xt_s2 .^2)+ (y0xt_s2 .^2))'; threadend; 
								threadbegin;allxthetas3 =              (y1xt3 - y0xt3)';    allxstds3 =              sqrt((y1xt_s3 .^2)+ (y0xt_s3 .^2))'; threadend;  
								threadjoin;
								threadbegin;allxthetas_m  =            (y1xt_m  - y0xt_m)';     allxstds_m  =              sqrt((y1xt_s_m .^2) + (y0xt_s_m  .^2))'; threadend; 
								threadbegin;allxthetas2_m =            (y1xt2_m - y0xt2_m)';    allxstds2_m =              sqrt((y1xt_s2_m .^2)+ (y0xt_s2_m .^2))'; threadend; 
								threadbegin;allxthetas3_m =            (y1xt3_m - y0xt3_m)';    allxstds3_m =              sqrt((y1xt_s3_m .^2)+ (y0xt_s3_m .^2))'; threadend;  
								threadjoin;
		
							else; 
								threadbegin;allxthetas  = allxthetas  | ((y1xt - y0xt)');   allxstds  = allxstds  | (sqrt((y1xt_s .^2) + (y0xt_s  .^2))'); threadend; 
								threadbegin;allxthetas2 = allxthetas2 | ((y1xt2 - y0xt2)'); allxstds2 = allxstds2 | (sqrt((y1xt_s2 .^2)+ (y0xt_s2 .^2))'); threadend; 
								threadbegin;allxthetas3 = allxthetas3 | ((y1xt3 - y0xt3)'); allxstds3 = allxstds3 | (sqrt((y1xt_s3 .^2)+ (y0xt_s3 .^2))'); threadend; 
								threadjoin;	
								threadbegin;allxthetas_m  = allxthetas_m  | ((y1xt_m - y0xt_m)');   allxstds_m  = allxstds_m  | (sqrt((y1xt_s_m .^2) + (y0xt_s_m  .^2))'); threadend; 
								threadbegin;allxthetas2_m = allxthetas2_m | ((y1xt2_m - y0xt2_m)'); allxstds2_m = allxstds2_m | (sqrt((y1xt_s2_m .^2)+ (y0xt_s2_m .^2))'); threadend; 
								threadbegin;allxthetas3_m = allxthetas3_m | ((y1xt3_m - y0xt3_m)'); allxstds3_m = allxstds3_m | (sqrt((y1xt_s3_m .^2)+ (y0xt_s3_m .^2))'); threadend; 
								threadjoin;			
			
							endif;
							if b_i == 0;
								" ========================  treatment ";; format 2,0; i;; format 8,8; $treatment;
								" ========================  population ";; format 2,0; ate_i;; format 8,8; $treatment;
								" ---------------------------  P(X) ATENT ------NONTREATED ----- P(X) -------------";
								"============== P(X) ATENT ======== NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
								call nice_printing_binary_new(outdat_t3,all_outcomes,allxthetas3,allxstds3); 
								"============== P(X) ATENT ======== Linear bias correction (weights used for inference take bias adjustment into account) ========";
								call nice_printing_binary_new(outdat_t1,all_outcomes,allxthetas,allxstds);
								"============== P(X) ATENT =======  This version uses weighted logits instead of weighted ols for the binary outcomes (weights used for inference take bias adjustment NOT into account)====================";
								call nice_printing_binary_new(outdat_t2,all_outcomes,allxthetas2,allxstds2);  
								" --------------------------- P(X,M) ATENT ------NONTREATED ----- P(X,M) -------------";
								"=============== P(X,M) ATENT ======= NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
								call nice_printing_binary_new(outdat_t3_M,all_outcomes,allxthetas3_M,allxstds3_M); 
								"================P(X,M) ATENT ====== Linear bias correction (weights used for inference take bias adjustment into account) ========";
								call nice_printing_binary_new(outdat_t1_M,all_outcomes,allxthetas_M,allxstds_M);
								"===============P(X,M) ATENT======  This version uses weighted logits instead of weighted ols for the binary outcomes (weights used for inference take bias adjustment NOT into account)====================";
								call nice_printing_binary_new(outdat_t2_M,all_outcomes,allxthetas2_M,allxstds2_M);  
		
							endif;
						else;   /* ------------------- ATET -------------------------- */
						outdat_t4 = outdat4 $+ "_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file f�r alle populationen */
						outdat_t5 = outdat5 $+ "_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file f�r alle populationen */
						outdat_t6 = outdat6 $+ "_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file f�r alle populationen */
						outdat_t4_M = outdat4 $+ "_M_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file f�r alle populationen */
						outdat_t5_M = outdat5 $+ "_M_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file f�r alle populationen */
						outdat_t6_M = outdat6 $+ "_M_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file f�r alle populationen */
						if no_w;
							threadstat {y0_t, y1_t, y0_t_s, y1_t_s, y1_t_975, y0_t_975, y1_t_025, y0_t_025, n0, n1, n_ref,w_match_WOLS}	=evaluate_m_s_all_2(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,1,t_val_k); 
							threadstat {y0_t2,y1_t2,y0_t_s2,y1_t_s2,y1_t_9752,y0_t_9752,y1_t_0252,y0_t_0252,n02,n12,n_ref2,w_match2}		=evaluate_m_s_all_2(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,2,t_val_k); 
							threadstat {y0_t3,y1_t3,y0_t_s3,y1_t_s3,y1_t_9753,y0_t_9753,y1_t_0253,y0_t_0253,n03,n13,n_ref3,w_match3}		=evaluate_m_s_all_2(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,0,t_val_k);  
							threadjoin;	
							threadstat {y0_t_M, y1_t_M, y0_t_s_M, y1_t_s_M, y1_t_975_M, y0_t_975_M, y1_t_025_M, y0_t_025_M, n0_M, n1_M, n_ref_M,w_match_WOLS_M}	=evaluate_m_s_all_2(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,1,t_val_k); 
							threadstat {y0_t2_M,y1_t2_M,y0_t_s2_M,y1_t_s2_M,y1_t_9752_M,y0_t_9752_M,y1_t_0252_M,y0_t_0252_M,n02_M,n12_M,n_ref2_M,w_match2_M}	=evaluate_m_s_all_2(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,2,t_val_k); 
							threadstat {y0_t3_M,y1_t3_M,y0_t_s3_M,y1_t_s3_M,y1_t_9753_M,y0_t_9753_M,y1_t_0253_M,y0_t_0253_M,n03_M,n13_M,n_ref3_M,w_match3_M}	=evaluate_m_s_all_2(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,0,t_val_k);  
							threadjoin;					
						else;	
							threadstat {y0_t, y1_t, y0_t_s, y1_t_s, y1_t_975, y0_t_975, y1_t_025, y0_t_025, n0, n1, n_ref,w_match_WOLS,w1}	=evaluate_m_s_all_2w(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,1,sample_w,t_val_k); 
							threadstat {y0_t2,y1_t2,y0_t_s2,y1_t_s2,y1_t_9752,y0_t_9752,y1_t_0252,y0_t_0252,n02,n12,n_ref2,w_match2,w12}	=evaluate_m_s_all_2w(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,2,sample_w,t_val_k); 
							threadstat {y0_t3,y1_t3,y0_t_s3,y1_t_s3,y1_t_9753,y0_t_9753,y1_t_0253,y0_t_0253,n03,n13,n_ref3,w_match3,w13}	=evaluate_m_s_all_2w(temp8,dep_name,all_outcomes,p_score_name_t,x_extra_name,w_name,0,sample_w,t_val_k);  
							threadjoin;	
							threadstat {y0_t_M, y1_t_M, y0_t_s_M, y1_t_s_M, y1_t_975_M, y0_t_975_M, y1_t_025_M, y0_t_025_M, n0_M, n1_M, n_ref_M,w_match_WOLS_M,w1_M}=evaluate_m_s_all_2w(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,1,sample_w,t_val_k); 
							threadstat {y0_t2_M,y1_t2_M,y0_t_s2_M,y1_t_s2_M,y1_t_9752_M,y0_t_9752_M,y1_t_0252_M,y0_t_0252_M,n02_M,n12_M,n_ref2_M,w_match2_M,w12_M}	=evaluate_m_s_all_2w(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,2,sample_w,t_val_k); 
							threadstat {y0_t3_M,y1_t3_M,y0_t_s3_M,y1_t_s3_M,y1_t_9753_M,y0_t_9753_M,y1_t_0253_M,y0_t_0253_M,n03_M,n13_M,n_ref3_M,w_match3_M,w13_M}	=evaluate_m_s_all_2w(temp8_M,dep_name,all_outcomes,p_score_name_t_M,x_extra_name_M,w_name_M,0,sample_w,t_val_k);  
							threadjoin;	
						
						endif;
						if b_i == 0; 
							?;output on;
							"========================================================================";
							"====== Balancing Test before matching ATET ====";
							if no_w;	?;"P(X)";	call test_quality_of_matched_s_m_s(temp8,mcv|score_name_all,dep_Name,1); 
										?;"P(X,M)";	call test_quality_of_matched_s_m_s(temp8_M,mcv_m|score_name_all_M,dep_Name,1);
							else;		w0x = readvar(temp8,sample_w);dxx = readvar(temp8,dep_name);
										?;"P(X)";	call test_quality_of_matched_s_m_s_w(temp8,mcv|score_name_all,dep_Name,selif(w0x,dxx .== 0),w1); 
										w0x_M = readvar(temp8_M,sample_w);dxx_M = readvar(temp8_M,dep_name);
										?;"P(X,M)";	call test_quality_of_matched_s_m_s_w(temp8_M,mcv_M|score_name_all_M,dep_Name,selif(w0x_M,dxx_M .== 0),w1_M); 
							endif;	
							"========================================================================";
							"====== Balancing Test with weights from radius matching ATET ====";
											
							if no_w;	?;"P(X)";	call test_quality_of_matched_s_m_s(temp8,mcv|score_name_all,dep_Name,w_match3);
										?;"P(X,M)";	call test_quality_of_matched_s_m_s(temp8_M,mcv_M|score_name_all_M,dep_Name,w_match3_M); 	
							else;		?;"P(X)";	call test_quality_of_matched_s_m_s_w(temp8,mcv|score_name_all,dep_Name,w_match3,w13);
										?;"P(X,M)";	call test_quality_of_matched_s_m_s_w(temp8_M,mcv_M|score_name_all_M,dep_Name,w_match3_M,w13_M);
							endif;
							"========================================================================";
							"====== Balancing Test with weights from linear bias adjustment ATET ====";
							if no_w;	?;"P(X)";	call test_quality_of_matched_s_m_s(temp8,mcv|score_name_all,dep_Name,w_match_WOLS); 
										?;"P(X,M)";	call test_quality_of_matched_s_m_s(temp8_M,mcv_M|score_name_all_M,dep_Name,w_match_WOLS_M);
							else;		?;"P(X)";	call test_quality_of_matched_s_m_s_w(temp8,mcv|score_name_all,dep_Name,w_match_WOLS,w1); 
										?;"P(X,M)";	call test_quality_of_matched_s_m_s_w(temp8_M,mcv_M|score_name_all_M,dep_Name,w_match_WOLS_M,w1_M); 
							endif;		
    
							call write_outcomes_in_file(outdat_t4,y0_t, y1_t, y0_t_s, y1_t_s, y1_t_975, y0_t_975, y1_t_025, y0_t_025, n0, n1, n_ref,share_remaining,all_outcomes,create_yes); 
							call write_outcomes_in_file(outdat_t5,y0_t2,y1_t2,y0_t_s2,y1_t_s2,y1_t_9752,y0_t_9752,y1_t_0252,y0_t_0252,n02,n12,n_ref2,share_remaining,all_outcomes,create_yes); 
							call write_outcomes_in_file(outdat_t6,y0_t3,y1_t3,y0_t_s3,y1_t_s3,y1_t_9753,y0_t_9753,y1_t_0253,y0_t_0253,n03,n13,n_ref3,share_remaining,all_outcomes,create_yes); 
							
							call write_outcomes_in_file(outdat_t4_M,y0_t_M, y1_t_M, y0_t_s_M, y1_t_s_M, y1_t_975_M, y0_t_975_M, y1_t_025_M, y0_t_025_M, n0_M, n1_M, n_ref_M,share_remaining,all_outcomes,create_yes); 
							call write_outcomes_in_file(outdat_t5_M,y0_t2_M,y1_t2_M,y0_t_s2_M,y1_t_s2_M,y1_t_9752_M,y0_t_9752_M,y1_t_0252_M,y0_t_0252_M,n02_M,n12_M,n_ref2_M,share_remaining,all_outcomes,create_yes); 
							call write_outcomes_in_file(outdat_t6_M,y0_t3_M,y1_t3_M,y0_t_s3_M,y1_t_s3_M,y1_t_9753_M,y0_t_9753_M,y1_t_0253_M,y0_t_0253_M,n03_M,n13_M,n_ref3_M,share_remaining,all_outcomes,create_yes); 
	
						endif;
						if (b_i == 0) or (b_i == 1); 
							threadbegin; all_thetas  = 	(y1_t - y0_t)';      all_stds  =	sqrt((y1_t_s .^2) + (y0_t_s  .^2))';	threadend;
							threadbegin; all_thetas2 =  (y1_t2 - y0_t2)';    all_stds2 =    sqrt((y1_t_s2 .^2)+ (y0_t_s2 .^2))';	threadend;
							threadbegin; all_thetas3 =  (y1_t3 - y0_t3)';    all_stds3 =    sqrt((y1_t_s3 .^2)+ (y0_t_s3 .^2))'; 	threadend; 
							threadjoin;	
							
							threadbegin; all_thetas_M  =(y1_t_M - y0_t_M)';  all_stds_M  =  sqrt((y1_t_s_M .^2) + (y0_t_s_M  .^2))';	threadend;
							threadbegin; all_thetas2_M =(y1_t2_M - y0_t2_M)';all_stds2_M =  sqrt((y1_t_s2_M .^2)+ (y0_t_s2_M .^2))';	threadend;
							threadbegin; all_thetas3_M =(y1_t3_M - y0_t3_M)';all_stds3_M =  sqrt((y1_t_s3_M .^2)+ (y0_t_s3_M .^2))'; 	threadend; 
							threadjoin;					
						else; 
							threadbegin;all_thetas  = all_thetas | ((y1_t - y0_t)');    all_stds  = all_stds  | (sqrt((y1_t_s .^2) + (y0_t_s  .^2))');	threadend;
							threadbegin;all_thetas2 = all_thetas2 |((y1_t2 - y0_t2)');  all_stds2 = all_stds2 | (sqrt((y1_t_s2 .^2)+ (y0_t_s2 .^2))');	threadend;
							threadbegin;all_thetas3 = all_thetas3 |((y1_t3 - y0_t3)');  all_stds3 = all_stds3 | (sqrt((y1_t_s3 .^2)+ (y0_t_s3 .^2))');	threadend;    
							threadjoin;	
							threadbegin;all_thetas_M  = all_thetas_M | ((y1_t_M - y0_t_M)');    all_stds_M  = all_stds_M  | (sqrt((y1_t_s_M .^2) + (y0_t_s_M  .^2))');	threadend;
							threadbegin;all_thetas2_M = all_thetas2_M |((y1_t2_M - y0_t2_M)');  all_stds2_M = all_stds2_M | (sqrt((y1_t_s2_M .^2)+ (y0_t_s2_M .^2))');	threadend;
							threadbegin;all_thetas3_M = all_thetas3_M |((y1_t3_M - y0_t3_M)');  all_stds3_M = all_stds3_M | (sqrt((y1_t_s3_M .^2)+ (y0_t_s3_M .^2))');	threadend;    
							threadjoin;	
						endif;
						
						if b_i == 0; 
								
							outdat_t7_M = outdat7 $+ "__M_ID_T_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */
							outdat_t8_M = outdat8 $+ "_M_ID_NT" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */

							indirect_T = all_thetas - all_thetas_M;  y1idt = y0xt + indirect_T'; 	y0idt = y0xt;
							indirect_NT = allxthetas - allxthetas_M; y0idNt = y1xt - indirect_NT';	y1idNt = y1xt;
							
							call write_outcomes_in_file(outdat_t7_M,y0IDt, y1IDt, y0_t_s, y1_t_s,y1_t_975,y0_t_975,y1_t_025,y0_t_025,n0,n1,n_ref,share_remaining,all_outcomes,create_yes); 
							call write_outcomes_in_file(outdat_t8_M,y0IDNt,y1IDNt,y0_t_s, y1_t_s,y1_t_975,y0_t_975,y1_t_025,y0_t_025,n0,n1,n_ref,share_remaining,all_outcomes,create_yes); 
						endif;
						
						
						
						if b_i == 0;
							?;" --------------------------------------------------------------------------------------------------------------------------";
							"========================  population ";; format 2,0; ate_i;; format 8,8; $treatment;
							" ------------P(X) --------------- P(X) ATET ------ TREATED -----  -------------";
							"================P(X) ATET====== NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
							call nice_printing_binary_new(outdat_t6,all_outcomes,all_thetas3,all_stds3); 
							"================P(X) ATET====== Linear bias correction (weights used for inference take bias adjustment into account) ========";
							call nice_printing_binary_new(outdat_t4,all_outcomes,all_thetas,all_stds);
							"================P(X) ATET======  Bias correction based on weighted logits instead of weighted ols for binary outcomes (weights used for inference take bias adjustment NOT into account)";
							call nice_printing_binary_new(outdat_t5,all_outcomes,all_thetas2,all_stds2);  
							
							?;" --------------------------------------------------------------------------------------------------------------------------";
							"========================  population ";; format 2,0; ate_i;; format 8,8; $treatment;
							" ------------P(X,M) --------------- P(X,M) ATET ------ TREATED -----  -------------";
							"================P(X,M) ATET====== NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
							call nice_printing_binary_new(outdat_t6,all_outcomes,all_thetas3_M,all_stds3_M); 
							"================P(X,M) ATET====== Linear bias correction (weights used for inference take bias adjustment into account) ========";
							call nice_printing_binary_new(outdat_t4,all_outcomes,all_thetas_M,all_stds_M);
							"================P(X,M) ATET======  Bias correction based on weighted logits instead of weighted ols for binary outcomes (weights used for inference take bias adjustment NOT into account)";
							call nice_printing_binary_new(outdat_t5,all_outcomes,all_thetas2_M,all_stds2_M);  
							
						endif;  
						save allxthetas_M,allxthetas2_M,allxthetas3_M,allxstds_M,allxstds2_M,allxstds3_M;
						save allxthetas,allxthetas2,allxthetas3,allxstds,allxstds2,allxstds3; 	
						save all_thetas_M,all_thetas2_M,all_thetas3_M,all_stds_M,all_stds2_M,all_stds3_M;
						save all_thetas,all_thetas2,all_thetas3,all_stds,all_stds2,all_stds3; 
	

/*
						/* --------------------------------- ATE ------------------------------------ */
						outdat_t7 = outdat7 $+ "_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */
						outdat_t8 = outdat8 $+ "_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */
						outdat_t9 = outdat9 $+ "_" $+ ftos(i,"%0*.*lf",2,0) $+ "_" $+ ftos(ate_i,"%0*.*lf",1,0); /* immer das gleiche output file fuer alle populationen */
						threadbegin;y0at  	= 			n0  ./ (n1 + n0)  	.* y0xt  		+ 	n1  ./ (n1 + n0)  	.* y0_t; 
									y1at  	= 			n0  ./ (n1 + n0)  	.* y1xt  		+ 	n1  ./ (n1 + n0)  	.* y1_t;
									y0at_s  = sqrt(  (	n0  ./ (n1 + n0)  	.* y0xt_s) .^2  + (	n1  ./ (n1 + n0)  	.* y0_t_s) .^2 );    
									y1at_s  = sqrt(  (	n0  ./ (n1 + n0)  	.* y1xt_s) .^2  + (	n1  ./ (n1 + n0)  	.* y1_t_s) .^2 );  
									
									y0at2 	= 			n02 ./ (n12+ n02) 	.* y0xt2		+ 	n12 ./ (n12+ n02) 	.* y0_t2;  	
									y1at2 	= 			n02 ./ (n12+ n02) 	.* y1xt2 		+ 	n12 ./ (n12+ n02) 	.* y1_t2;
						threadend;				
						threadbegin;y0at_s2 = sqrt(  (	n02  ./ (n12 + n02) .* y0xt_s2).^2  + (	n12 ./ (n12 + n02) 	.* y0_t_s2).^2 );    
									y1at_s2 = sqrt(  (	n02  ./ (n12 + n02) .* y1xt_s2).^2  + (	n12 ./ (n12 + n02) 	.* y1_t_s2).^2 );					
									
									y0at3 	= 			n03 ./ (n13+ n03) 	.* y0xt3 		+ 	n13 ./ (n13+ n03) 	.* y0_t3; 
									y1at3 	= 			n03 ./ (n13+ n03) 	.* y1xt3 		+ 	n13 ./ (n13+ n03) 	.* y1_t3;             
									y0at_s3 = sqrt(  (	n03  ./ (n13 + n03) .* y0xt_s3).^2  + (	n13 ./ (n13 + n03) 	.* y0_t_s3).^2 ); 
									y1at_s3 = sqrt(  (	n03  ./ (n13 + n03) .* y1xt_s3).^2  + (	n13 ./ (n13 + n03) 	.* y1_t_s3).^2 );
						threadend;			
						threadjoin;	
						if b_i == 0; 
							call write_outcomes_in_file(outdat_t7,y0at, y1at, y0at_s, y1at_s, y1_t_975, y0_t_975, y1_t_025, y0_t_025, n0, n1, n_ref,share_remaining,all_outcomes,create_yes); 
							call write_outcomes_in_file(outdat_t8,y0at2,y1at2,y0at_s2,y1at_s2,y1_t_9752,y0_t_9752,y1_t_0252,y0_t_0252,n02,n12,n_ref2,share_remaining,all_outcomes,create_yes); 
							call write_outcomes_in_file(outdat_t9,y0at3,y1at3,y0at_s3,y1at_s3,y1_t_9753,y0_t_9753,y1_t_0253,y0_t_0253,n03,n13,n_ref3,share_remaining,all_outcomes,create_yes); 
						endif;
						if (b_i == 0) or (b_i == 1); 
							threadbegin; allathetas  =              (y1at - y0at)';      allastds  =              sqrt((y1at_s .^2) + (y0at_s  .^2))';threadend;
							threadbegin; allathetas2 =              (y1at2 - y0at2)';    allastds2 =              sqrt((y1at_s2 .^2)+ (y0at_s2 .^2))';threadend;  
							threadbegin; allathetas3 =              (y1at3 - y0at3)';    allastds3 =              sqrt((y1at_s3 .^2)+ (y0at_s3 .^2))';threadend; 
							threadjoin;	
						else;                                                                                                            
							threadbegin; allathetas  = allathetas  | ((y1at - y0at)');   allastds  = allastds  | (sqrt((y1at_s .^2) + (y0at_s  .^2))');threadend;
							threadbegin; allathetas2 = allathetas2 | ((y1at2 - y0at2)'); allastds2 = allastds2 | (sqrt((y1at_s2 .^2)+ (y0at_s2 .^2))');threadend;
							threadbegin; allathetas3 = allathetas3 | ((y1at3 - y0at3)'); allastds3 = allastds3 | (sqrt((y1at_s3 .^2)+ (y0at_s3 .^2))');threadend;
							threadjoin;	
						endif;
						
						if b_i == 0;
							?;" --------------------------------------------------------------------------------------------------------------------------";
							" ========================  population ";; "all";; format 8,8; $treatment;
							" ---------------------------  ATE ------ population ------------------";
							"====================== NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
							call nice_printing_binary_new(outdat_t9,all_outcomes,all_thetas3,all_stds3); 
							"====================== Linear bias correction (weights used for inference take bias adjustment into account) ========";
							call nice_printing_binary_new(outdat_t7,all_outcomes,allathetas,allastds);
							"=====================  This version uses weighted logits instead of weighted ols for the binary outcomes (weights used for inference take bias adjustment NOT into account)====================";
							call nice_printing_binary_new(outdat_t8,all_outcomes,allathetas2,allastds2);       
						endif;
ATE */						
					endif; 
				endfor; /* ate_i */
			endif; @ i /= j @
			endfor; /* j */
		endfor; /* i */
	endif;
endfor; @ bootstrap @



if bootstrap > 1; output on;
  ?;"FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFinal resultsssssssssssssssssssssssssssssssssssssssssssss"; 
  "Common support and excess weight adjustments are based on p(x,m) [not on p(x)]";	
  "---------------------------------------------------------------------------------";
  "-------------------------------- P(X) ATENT  NONTREATED -----------------------------------";
  "====================P(X) ATENT == NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
  call nice_printing_binary_new_2(outdat_t3,all_outcomes,allxthetas3,allxstds3,p_val_kernel); 
  "====================P(X) ATENT== Linear bias correction (weights used for inference take bias adjustment into account) ========";
  call nice_printing_binary_new_2(outdat_t1,all_outcomes,allxthetas,allxstds,p_val_kernel);
  "=  This version uses weighted logits instead of weighted ols for the binary outcomes (weights used for inference take bias adjustment NOT into account) =";
  call nice_printing_binary_new_2(outdat_t2,all_outcomes,allxthetas2,allxstds2,p_val_kernel);  
  ?;?;"---------------------------------------------------------------------------------";
  "--------------------------------- P(X) ATET  TREATED -----------------------------------";
  "====================== P(X) ATET NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
  call nice_printing_binary_new_2(outdat_t6,all_outcomes,all_thetas3,all_stds3,p_val_kernel); 
  "====================== P(X) ATET Linear bias correction (weights used for inference take bias adjustment into account) ========";
  call nice_printing_binary_new_2(outdat_t4,all_outcomes,all_thetas,all_stds,p_val_kernel);
  "=  P(X) ATET Weighted logits instead of weighted ols for the binary outcomes (weights used for inference take bias adjustment NOT into account)=";
  call nice_printing_binary_new_2(outdat_t5,all_outcomes,all_thetas2,all_stds2,p_val_kernel);  
  ?;?;"---------------------------------------------------------------------------------";
  
  "-------------------------------- P(X,M) ATENT  NONTREATED -----------------------------------";
  "====================P(X,M) ATENT == NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
  call nice_printing_binary_new_2(outdat_t3_M,all_outcomes,allxthetas3_M,allxstds3_M,p_val_kernel); 
  "====================P(X,M) ATENT== Linear bias correction (weights used for inference take bias adjustment into account) ========";
  call nice_printing_binary_new_2(outdat_t1_M,all_outcomes,allxthetas_M,allxstds_M,p_val_kernel);
  "=  This version uses weighted logits instead of weighted ols for the binary outcomes (weights used for inference take bias adjustment NOT into account) =";
  call nice_printing_binary_new_2(outdat_t2_M,all_outcomes,allxthetas2_M,allxstds2_M,p_val_kernel);  
  ?;?;"---------------------------------------------------------------------------------";
  "--------------------------------- P(X,M) ATET  TREATED -----------------------------------";
  "====================== P(X,M) ATET NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
  call nice_printing_binary_new_2(outdat_t6_M,all_outcomes,all_thetas3_M,all_stds3_M,p_val_kernel); 
  "====================== P(X,M) ATET Linear bias correction (weights used for inference take bias adjustment into account) ========";
  call nice_printing_binary_new_2(outdat_t4_M,all_outcomes,all_thetas_M,all_stds_M,p_val_kernel);
  "=  P(X) ATET Weighted logits instead of weighted ols for the binary outcomes (weights used for inference take bias adjustment NOT into account)=";
  call nice_printing_binary_new_2(outdat_t5_M,all_outcomes,all_thetas2_M,all_stds2_M,p_val_kernel);  
  ?;?;"---------------------------------------------------------------------------------";
  /*
  " ========================  population ";; "all";; format 8,8; $treatment;
  " ---------------------------  ATE ------ population ------------------";
  "Asymptotic standard errors computed from weighted average of variances of ATET and ATENT; i.e. "
  "they are computed under the simplifying assumption that Cov(ATET, ATENT) = 0; "
  "Compare bootstrap variances to asymptotic variance to check quality of this approximation";
  "====================== NO bias correction (weights used for inference take bias adjustment NOT into account)========"; 
  call nice_printing_binary_new_2(outdat_t9,all_outcomes,allathetas3,allastds3,p_val_kernel); 
  "====================== Linear bias correction (weights used for inference take bias adjustment into account) ========";
  call nice_printing_binary_new_2(outdat_t7,all_outcomes,allathetas,allastds,p_val_kernel);
  "=  This version uses weighted logits instead of weighted ols for the binary outcomes (weights used for inference take bias adjustment NOT into account)=";
  call nice_printing_binary_new_2(outdat_t8,all_outcomes,allathetas2,allastds2,p_val_kernel);       
	*/
  ?;	
  "MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM";	
  "Mediation analysis	- all effects based on linear bias adjustment (check results above for robustness) ";
  "MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM";
  
  "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT";
  "TOTAL effect on the treated E[E(Y1|X,D=1)-E(Y0|X,D=1)] = ATET(x)";
   "------------------------------------------------------------------";
  "====================== P(X) ATET Linear bias correction (weights used for inference take bias adjustment into account) ========";
  call nice_printing_binary_new_2(outdat_t4,all_outcomes,all_thetas,all_stds,p_val_kernel);
  "-------------------------------------------------------------------";
  "DIRECT effect of the treatment on the treated under treatment E[E(Y1|X,M,D=1)-E(Y0|X,M,D=1)] = ATET(X,M)";
  "-------------------------------------------------------------------";
  "====================== P(X,M) ATET Linear bias correction (weights used for inference take bias adjustment into account) ========";
  call nice_printing_binary_new_2(outdat_t4_M,all_outcomes,all_thetas_M,all_stds_M,p_val_kernel);
  "-------------------------------------------------------------------";
  "INDIRECT effect of the treatment on the treated under no treatment ATET(X) - ATET(X,M)";
  "-------------------------------------------------------------------";
  "P-values based on bootstrapped t-statistic computed assuming ATET(X) and ATET(X,M) are independent";
  indirect_T = all_thetas - all_thetas_M; std_indirect_T = sqrt((all_stds .^2) + (all_stds_M .^2));
  call nice_printing_binary_new_2(outdat_t7_M,all_outcomes,indirect_T,std_indirect_T,p_val_kernel);
  "---------------------------------------------------------------------------------------------------------------------";
  
  "NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN";
  "TOTAL effect on the nontreated E[E(Y1|X,D=0)-E(Y0|X,D=0)] = ATENT(x)";
  "------------------------------------------------------------------";
  "====================P(X) ATENT== Linear bias correction (weights used for inference take bias adjustment into account) ========";
  call nice_printing_binary_new_2(outdat_t1,all_outcomes,allxthetas,allxstds,p_val_kernel);
  "-------------------------------------------------------------------";
  "DIRECT effect of the treatment on the untreated under no treatment E[E(Y1|X,M,D=0)-E(Y0|X,M,D=0)] = ATENT(X,M)";
  "-------------------------------------------------------------------";
  "====================P(X,M) ATENT== Linear bias correction (weights used for inference take bias adjustment into account) ========";
  call nice_printing_binary_new_2(outdat_t1_M,all_outcomes,allxthetas_M,allxstds_M,p_val_kernel);
   "-------------------------------------------------------------------";
  "INDIRECT effect of the treatment on the untreated under treatment ATENT(X) - ATENT(X,M)";
  "P-values based on bootstrapped t-statistic computed assuming ATENT(X) and ATENT(X,M) are independent";
  "-------------------------------------------------------------------";
  indirect_NT = allxthetas - allxthetas_M; std_indirect_NT = sqrt((allxstds .^2) + (allxstds_M .^2));
  call nice_printing_binary_new_2(outdat_t8_M,all_outcomes,indirect_NT,std_indirect_NT,p_val_kernel);
  
endif;
output on; call prog_end(date~time,ti1);
end;


@ ========================================================================================================= @
@ ========================================================================================================= @
/*
#include c:\mlechner\mediation\inc\probit_logit.inc;  
#include c:\mlechner\mediation\inc\matching.inc;  
#include c:\mlechner\mediation\inc\common_support.inc;  
#include c:\mlechner\mediation\inc\statistic_tools.inc;
#include c:\mlechner\mediation\inc\general_purpose.inc;
#include c:\mlechner\mediation\inc\bootstrap.inc;
*/
/**/
#include d:\mlechner\mediation\inc\probit_logit.inc;  
#include d:\mlechner\mediation\inc\matching.inc;  
#include d:\mlechner\mediation\inc\common_support.inc;  
#include d:\mlechner\mediation\inc\statistic_tools.inc;
#include d:\mlechner\mediation\inc\general_purpose.inc;
#include d:\mlechner\mediation\inc\bootstrap.inc;
/**/
@ ========================================================================================================= @
@ ========================================================================================================= @ 



